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We generalize the treatment of the electronic spin degrees of freedom in density functional calcu- 
lations to the case where the spin vector variables employed in the definition of the energy functional 
can vary in any direction in space. The expression for the generalized exchange-correlation potential 
matrix elements is derived for general functionals which among their ingredients include the electron 
density, its gradient and Laplacian, the kinetic energy density, and non-local Hartree-Fock type ex- 
change. We present calculations on planar Cr clusters that exhibit ground states with noncollinear 
spin densities due to geometrically frustrated antiferromagnetic interactions. 
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I. INTRODUCTION 

Since the spin polarized formulation 1 of density func- 
tional theory (DFT) 2 -^ was introduced, a large number of 
applications have been carried out in magnetic systems. 
In most cases, the spin density is assumed to adopt a 
single direction (collinear) at each point in space that is 
usually taken as z. However, there are a number of sys- 
tems where the spin density (or magnetization density) 
can take a more complicated structure, and vary its di- 
rection at each point in space. These type of noncollinear 
structures were observed in the form of helical spin den- 
sity waves or spin spirals for the ground state of 7-Fe^^ 
in geometrically frustrated systems like for instance the 
Kagome antiferromagnetic lattice^, and in systems with 
competing magnetic interactions such as the composite 
magnet LaMn 2 Gez£ and Fe . 5 Co . 5 Si4 

Several papers dealing with noncollinear spin density 
in DFT calculations have been published in the litera- 
ture. The pioneer work of Kubler and coworkers^ for the 
noncollinear local spin density approximation (LSDA) 
was later followed by a number of independent imple- 
mentations and applications. Most of these implemen- 
tations were carried out using periodic boundary con- 
ditions and plane waves, and are based either on the 
LSDA 10,11 or on a generalized gradient approximation 
(GGA). 12 Yamanaka and coworkers developed a general- 
ized DFT code based on Gaussian type orbitals.— Some 
noncollinear DFT calculations have been published deal- 
ing with magnetic crystals and with fourth period 
transition metal clusters i 11 ' 14 ! 15 ! 16 In all cases, the real- 
ization of the LSDA and GGA employed in noncollinear 
calculations is the same as that developed for collinear 
spin systems. 

Different parametrizations of the exchange-correlation 
energy (E xc ) have been proposed beyond the LSDA and 
the GGA, incorporating more ingredients in the defini- 
tion of E xc . The third rung in this hierarchy^ is the 
meta-GGA, which includes the kinetic energy density as 
a functional ingredient. Also, hybrid density functionals, 
which contain a portion of Hartree-Fock type exchange, 



can be regarded as belonging to the fourth rung (hyper- 
GGA) in this pictured The purpose of this paper is to 
provide a consistent generalization for the treatment of 
noncollinear spin variables in DFT calculations beyond 
the LSDA. 



II. THEORY 

To allow for noncollinear spin states in density 
functional calculations, we start by introducing two- 
component spinors as Kohn-Sham (KS) orbitals: 



(1) 



where ipf and tpf are spatial orbitals that can be ex- 
panded in a linear combination of atomic orbitals, 

^f(r)=5>^(r) (<r = <*,0). (2) 



Using the KS formulation, the electronic energy is par- 
titioned into four contributions: 



E — Ei 



E 



N 



E.j + E x 



(3) 



where Et is the kinetic energy, En is the nuclear-electron 
interaction energy, E.j is the classical electron-electron 
Coulomb repulsion energy, and E xc is the exchange- 
correlation (XC) energy. Searching for stationary so- 
lutions of E is equivalent to solving the KS equations, 
which in terms of two-component spinors ^ read: 



(T + V N + J + V xc )*i = e^i 



(4) 



where T = —1/2 V 2 is the kinetic energy operator, Vn is 
the external electron-nuclear potential, J is the Coulomb 
operator, and V xc is the exchange-correlation (XC) po- 
tential. We shall here refer to the KS equations in a 
two-component spinor basis as the generalized KS (GKS) 
equations. Since T, Vn, and J are diagonal in the two- 
dimensional spin space, the only term in Eq. [4] that cou- 
ples ipf and ipf? is V xc (the spin-orbit operator, present 
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in a relativistic Hamiltonian, also couples the two spinor 
components). The potential V xc depends on the choice 
of E xc and therefore the coupling between ipf and ipf in 
the non-relativistic GKS equations depends exclusively 
on E xc . 

Let us first recall the standard formulation of E xc com- 
monly employed in (collinear) unrestricted KS (UKS) cal- 
culations. The general expression of E xc for a hybrid case 
can be cast as:— 



E xc — aE°' 



El 



(1 - a)E* x 



(5) 



where E° FA and E° FA are the exchange-correlation con- 
tributions to the energy at some (semi)local density func- 
tional approximation (DFA), respectively, E" F is the 
Hartree-Fock type exchange energy, and a is a mixing pa- 
rameter (0 < a < 1). The functional forms of E X FA and 
E° FA (as well as the parameter a) depend, of course, on 
the choice of the functional employed in the actual calcu- 
lation. A general expression for E° FA = aE° FA + E° FA 
can be written as 



EE 



J d 3 rf(Q), 



(6) 



where Q is a set of variables (included in the definition 
of E xc ): 

Q = |n Q , np, Vn a , Vrip, r a , Tp, V 2 n Q V 2 n^| . (7) 

Here n a and np are the a and (3 electron densities, and r Q 
and Tp are the a and j3 kinetic energy densities, respec- 
tively, representing the "up" and "down" components 
along the z axis. 

On the other hand, in the noncollincar case, where the 
vector component of the local variables employed in the 
definition of E xc can point in any direction, E X FA can be 
generalized as follows, 



*&° = d A r f{Q) , 



(8) 



where 



Q = 



-}• (9) 



Here the subindices + and — refer to variables expressed 
in a local reference frame along the local spin quantiza- 
tion axis. The definition of these variables is given in 
detail in Section II. Note that by replacing Q by Q, we 
have only added degrees of freedom to the local variables 
in such a way that they are compatible with any arbitrary 
choice of the local spin axis. The dependence of / (and 
therefore E™°) on these variables remains unchanged. 

Two conditions must be satisfied by the set of variables 
in Q. First, we should recover the standard collinear case 
if we allow spin polarization only in one direction, and 
therefore replace the labels + and — in -E" c c by a and 
/?, respectively. In other words, noncollinear GKS so- 
lutions should coincide with collinear solutions obtained 



with the standard UKS approximation in cases where 
the ground state solution is collinear. Second, for spin- 
independent Hamiltonians, any arbitrary choice (other 
than z) of the spin quantization axis should leave the 
energy unchanged. This means that any rigid rotation 
of all local reference frames should not change the total 
energy. 

We therefore assume that E x ° depends on the local 
variables + and — in the same manner as in the standard 
collinear case, and that E™° must be invariant under 
rigid rotations of the spin quantization axis. This is, of 
course, not the most general possibility to define energy 
functionals for noncollinear magnetic systems^ 

For practical applications, it is necessary to evaluate 
the XC potential matrix to be employed in the solution 
of the GKS equations. These matrix elements in a set of 
localized orbitals {4>^} can be written as: 



(V, 



NC\ 

XC lliv 



3„ df(Q) 



dP^u 
3 df(Q)dq dp 



p,qeQ 



dq dp dPuv ' 



(10) 



where P^„ are matrix elements of the generalized density 
matrix, 



P„ 




(11) 



i£occ 



and p represents variables that are linear in P^ v . The 
derivative df /dq remains the same as in the collinear 
case. The rest of thisjiection is devoted to the definition 
of the variables q € Q in the local reference frame, and 
to the evaluation of (V^ C ) M „. 



A. Density 

The main ingredient for the construction of E x ° in the 
LSDA is the generalized density, n, which can be written 
in a two-component spin space as 



- 1/ x 1 

n= -(n + m-a- ) = - 



n + m z 



in. 



n — m z 



,(12) 



where n is the electron density, m = (m Xl m yi m z ) 
is the spin density (magnetization) vector, and er = 
(p x , <j y , a z ) are the Pauli matrices. In terms of two- 
component spinors ^i, n and m are defined as 



n(r) = J2 *!( r )*»( r ) 

i£occ 



and 



m k (r) = 



J2 ^(rH^r) 



(k = x,y, 



(13) 



(14) 
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Using Eqs.[T][2l and fTTl it is straightforward to express 
n and m as a linear combination of matrix elements of 
the generalized density matrix P^ u . 

A local reference system where the generalized density 
n is diagonal can be obtained by rotating n into ft', 



n+ 
n- 



(15) 



where 



n ± = i(n ± m) = ^ ( n ± to 2 + to 2 + to 2 ) (16) 

are the eigenvalues of n. The densities n + and n_ can 
be regarded as the local analogous of n a and np. The 
potential V* c c is 



V 



SE NC 



1 



(f + + f-m-a 



(17) 



The matrix elements of V" c c are 



E / 



<9f 

Vj(<^„)d 3 r 



]=x,y,z 



E 

- E 

k—x,y,z 



of 



d(Vjm k ) 



<9m fc 



<Tk(<t>Md 3 r, (23) 



where the first term on the right-hand side of Eq. [23] 
contributes to the scalar potential £™ c °, and the second 
and last terms add to the magnetic potential B™°. The 
last term on the right-hand side of Eq. [53] arises from the 
fact that Vn± depends on m (Eq. [21]). Applying the 
chain rule, and making use of Eq. 1211 the derivatives of 
/ in Eq. [23] can be expressed as (we do not consider here 
derivatives of / arising from df /dn± since they where 
considered in the previous section) 



where 



f ± = 



df , df 



± 



dn^ dn_ ' 



(18) 



and m = hi/to is the unit vector in the direction of m. 
The matrix elements of V" c c can be evaluated straight- 
forwardly as 



(V£ a ) MV = / d 3 r^Vi 



(19) 



The potential V" c c in Eq.[T7]can be split into two con- 
tributions, 



V 



(20) 



where £ = f + /2 can be interpreted as a scalar (electro- 
static) potential and £?™ c c = /~m/2 as a spin-dependent 
(magnetic) potential, which for the LSDA is always par- 
allel to m. It is worth mentioning that in the limit of no 
magnetization (to —> 0), /~ — > and therefore £>™ c c — > 0, 
recovering the non-magnetic case. 



B. Density gradients 

Let us consider next GGA energy functionals. In this 
case the new ingredients in E™° are Vn + and Vn_, 
whose Cartesian component j (j = x,y, z) is 



\7jn± = 



dn± 

= 2\ iU ± m ^ TOfcV i TO ' c ) ' ( 21 ) 

k— x,y ,z 

We note in passing that this family of density functionals 
usually depends on the gradient of the density through 
the auxiliary quantities 20 



Jab = Vn • Vn 6 a, b 



(22) 



df 



1 



df 1 TO fe _ 



9(Vj*mfe) 2m 3 



(24) 



(25) 



and 



df 



L—x,y,z 



1 fV/m fc rn fe 



9mfc , * — ' Z I m " to 



9t - Vim} , (26) 

TO^ J 



where we have defined 



3 fc = — r ± 



d(Vkn + ) 5(Vfcn_) 



(27) 



Alternatively, one can calculate V" c c as the functional 
derivative 



5E NC 

xc e- 

on 

5E™° 



6n 



ySE^l 



a,. 



(28) 



Using the chain rule for functional derivatives, SE"^ /Srrii 
can be written as 



8E™ = 1 8E™ 
Srrii 2 Sn + 



se: 



Sri- 



-)- 



(29) 



Multiplying Eq. [2H by integrating over all space 

and applying integration by parts, we obtain Eq. 1231 
Therefore, from our definition of -B" c c , the contribution 
to the XC magnetic field from Eq. [28] is always par- 
allel to m. However, this is not necessarily the case for a 
general form of a GGA functional, as shown by Capelle 
and coworkers. 21 

One issue that is worth addressing is how our formu- 
lation differs from previous noncollinear generalizations 
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of the GGA. In our case, we employ Vn ± as denned 
in Eq. [3TJ and therefore the ingredients used in E™£ 
are strictly the gradients of the quantities n ± . Other 
implementation a 12 ' 22 employ either Vto or the z com- an d 
ponent of the projection of Vm^ onto m, therefore im- 
posing the constraint of £>™ c c being parallel to m. Here 
as in the LSDA case, in the limit of no magnetization 
(to — > and Vto — > 0), — > for k = x,y,z and hence 
£5™ c c — > 0, recovering the non-magnetic case. where 



C. Kinetic energy density 



df_ 

du k 



mi, 
— —/ 

2m 



1 



df 



dmi. 2m 



-h 



E 

j=x,y,z 



m m k 

2 / ' 

TO « 



9f_ ± 9f_ 
dr. dr_ 



(36) 



(37) 



(38) 



To deal with kinetic energy density contributions, we 
can proceed in analogy to Sec. Ill Al and define a general- 
ized kinetic energy density: 



r = — (r + u • a 



t + u z u x — lUy 
2 V u x + iu v t - u z 



(30) 



where r and u can be written in terms of two-component 
spinors as 



-^ = 1 E (V^r^-V^Ci 



(31) 



and 

Uk 



W = \ E ( Vv M r )) tff ^ • ( fc = y> z * 32 ) 



Comparing f (Eq. [30]) and n (Eq. [T^J) one is tempted 
to define r+ and r- as the local eigenvalues of f. How- 
ever, this choice would lead to a different local reference 
frame than the one used for n, and therefore collinear 
solutions obtained in this way will not necessarily be the 
same as those obtained with standard unrestricted KS 
calculations. To avoid this problem, we have chosen the 
following definition for r+ and r-, 



r ± = i(r±m-u). 



(33) 



which is equivalent to locally projecting u onto the axis 
defined by m. Using this choice for r+ and r-, the con- 
tribution to the XC potential matrix elements can be 
written as 

(KrW - E / %(V 3 ^V^)d 3 r 

df 



J=x,y,z 



E 

j : k—x,y,z 

- E 

k—x,y,z 



du k 
df 



o'fe(V ; ,-0 M Vj^ !/ )(i ,3 r 



a k (^)d 3 r. (34) 



The partial derivatives of / in Eq. [34] can be expressed 
in terms of the derivatives of / with respect to r+ and 
r- as 



dr 2 ' 



(35) 



The first term on the right-hand side of Eq. [M] con- 
tributes to the scalar potential £™ c c , while the second and 
third terms are spin dependent and therefore contribute 
to £>™ c c . In the limit where to — » then |m ■ u| — > 0, and 



therefore 23™ c c — > since ft- ->0. 



D. Laplacian of the density 

The dependence of E°^ A with the Laplacian of the 
density can be generalized for the noncollinear case 
through V 2 n±, which in terms of n, m, and their deriva- 
tives can be written as: 

W = ilv 2 n± y r^v 2 TO fc+ y (^^ 

2 I / — ' TO z — ' m 



k—x,y^z j—x,y,z 



l=x,y,z 



nr 



■■))}■ 



(39) 



Three types of contributions to the XC potential arise 
in this case, since V 2 n + and V 2 n_ depend on V 2 n, 
V 2 TOfe, mfc, and Vrrifc (Eq. 



(K™ C W = J(t + +t-a • rh)V 2 (<^)d 3 r 



]=x,y,z 



E 



rcr fc <9(V 2 m) 



k—x,y,z 



2 d(V jtiik 
t-o k d{V 2 m) 



k,j—x,y,z ' 

+ E /^^^V->""< <40) 



2 9TOfc 



where 



df ± df 



d(\7 2 n+) <9(V 2 n-) 



(41) 



The derivatives of V 2 to with respect to the linear vari- 
ables can be obtained from Eq. [3^1 In the limit of no 
magnetization where to — > and V 2 to — > 0, all contribu- 
tions to £>™ c ° are zero since i~ — > 0. 
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E. Hartree-Fock type exchange 

The HF type exchange contribution to {V xc )^ Vl needed 
for hybrid density functional calculations, can be gener- 
alized in terms of two-component spinors as 

^' = E p aY'(/*)> (42) 

A? 

where are the spin blocks of the generalized den- 

sity matrix in Eq. (|lip . The notation (/xi/|£A) has been 
introduced for the two electron integrals in the AO ba- 
sis set. This is analogous to the generalized unrestricted 
HF (GUHF) approximation^ The evaluation of the four 
matrix blocks of K™ is carried out by splitting P££ 
into real and imaginary parts, and symmetric and anti- 
symmetric components. The symmetric imaginary and 
antisymmetric real contributions to K aa and are 
zero because of the hermiticity requirement of the Kohn- 
Sham (or Hartree-Fock) Hamiltonian. Therefore, a total 
of eight HF exchange blocks need to be computed, four 
of them are symmetric and four are antisymmetric. 

III. IMPLEMENTATION 

We have implemented the SCF solution of the GKS 
equations in the Gaussian suite of programs^* Molecular 
spinors (Eq. [2]) are spanned in terms of atomic Gaussian 
orbitals using a set of complex coefficients cf ^ . These co- 
efficients are employed to construct the generalized den- 
sity matrix P M „ (Eq. [TTj) . from which the Hartree-Fock 
type exchange matrix can be calculated, as well as all the 
variables needed for the numerical quadrature employed 
in the evaluation of {V xc c )^. To accelerate the SCF con- 
vergence, we have generalized the direct inversion of the 
iterative subspace (DIIS)^ and the energy-based DIIS 26 
techniques for two-component complex spinors. 

In a GKS calculation, the spin density of the system 
is fully unconstrained, and is thus allowed to change in 
any arbitrary spatial direction. For instance, in a sim- 
ple calculation of the (nonrelativistic) ground state of 
the hydrogen atom using LSDA, one can obtain an in- 
finite manifold of solutions with the same total energy 
but different orientation of the spin density. These solu- 
tions are just linear combinations of the two degenerate 
linearly independent solutions. 

We have verified for a representative sample of func- 
tionals that for cases with ground state collinear solu- 
tions, the GKS solution for different choices of the quan- 
tization axis gives the same total energy as in a collinear 
UKS calculation. We have also verified that the calcu- 
lated electric dipole moment evaluated as finite differ- 
ences agrees with the expectation value of the electric 
dipole operator, satisfying the Hcllmann-Feynman theo- 
rem. 

The third term on the right-hand side of Eq. [53] leads 
to instabilities in the numerical integration due to the 




FIG. 1: Schematic representation of the Cr3, Crs, Cr7, and 
Cri2 clusters employed in our tests. The arrows represent 
the magnetization orientation on each atom as qualitatively 
obtained in our calculations. For Cr3 and Crs, two different 
chiralities were considered. 

presence of V(mj/m), which is exactly zero for collinear 
spin densities but may present significant oscillations for 
spin densities that are slightly noncollinear. For some 
cases, this prevented converging the total SCF energy 
better than 10 -6 hartree, which is not enough for our 
standard accurate convergence criteria. To avoid this 
problem, we discard contributions to Eq. [23] from grid 
points where the magnetization helicity, defined as mh — 
m • (V x m), is less than a certain threshold (rrih = in 
collinear cases). In our tests, a cutoff value of mt <10 
worked reasonably well. 



IV. RESULTS 

In order to test our GKS code, we have chosen a set of 
planar Cr clusters where the ground state is expected 
to exhibit noncollinear spin density arising from geo- 
metrically frustrated antiferromagnetic coupling (as in a 
Heisenberg spin Hamiltonian model) between neighbor- 
ing Cr atoms. In Fig. [T] we show a scheme of the Cr3 
(C 3v ), Cr 5 (C 5v ) , Cr 7 (C 6ll ), and Cr i2 (C 6u ) clusters 
and their resulting magnetic structures obtained in this 
work. In all cases, we have set the Cr-Cr bond length to 
3.70 Bohr in our calculations. This allows us to compare 
the direct effect of each functional on the magnetization. 

All calculations were carried out using a Ne core 
energy-consistent relativistic effective core potential 
(RECP) from the Stuttgart /Cologne group (Ref. 
We have employed a polarized triple-^ Gaussian ba- 
sis set consisting of 8s7p6dlf functions contracted to 
6s5p3dl/^l Even though our implementation offers the 
possibility of including the spin-orbit operator, either us- 
ing RECPs or in an all-electron framework, we have cho- 
sen not to include the spin-orbit interaction in the present 
test calculations. Atomic magnetic moments are calcu- 
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TABLE I: Atomic magnetic moments (in Bohr magnetons, 
Hb) and (S 2 ) (in fi%) of Cr clusters calculated using different 
energy functionals. See Fig. [T]for a scheme of the spin density 
configurations. 



Cluster 


Property 






Method 






LSDA 


PBE 


TPSS 


PBEh 


GUHF 


Cr 3 (C 3 „) 


m 


1.44 


1.66 


1.93 


2.40 


2.95 




(S 2 ) 


3.25 


3.87 


4.71 


6.32 


8.11 


Cr 5 (C 5v ) 


m 


1.61 


1.84 


2.07 


2.47 


2.91 




(S 2 ) 


5.21 


6.21 


7.38 


9.78 


12.47 


Cr 7 (C 6 „) 


m c 


0.18 


0.21 


0.31 


2.09 


2.36 




m e 


0.24 


0.89 


1.29 


2.33 


2.87 




(deg.) 


143 


103 


100 


105 


107 




(S 2 ) 


2.08 


4.02 


5.80 


14.60 


22.57 


Cr 12 (C 6v ) 


rrii 


0.72 


0.84 


1.03 


1.86 


2.15 




m e 


1.24 


1.50 


1.72 


2.28 


2.80 




(S 2 ) 


5.73 


7.60 


9.67 


17.71 


22.82 



lated according to Mulliken population analysis. The ex- 
pectation value (S 2 ) is evaluated in all DFT cases as for 
the GUHF determinant. 

We have chosen one representative functional from 
each class of functionals discussed in Section II. For 
the LSDA, we employ LDA (Dirac) exchange and the 
parametrization of Wosko, Wilk, and Nusair— for cor- 
relation (SVWN5); for the GGA we use the functional 
of Perdew, Burke, and Ernzerhof (PBE); 2 - for the meta- 
GGA we use the functional developed by Tao, Perdew, 
Staroverov, and Scuseria (TPSS)j^ and as a representa- 
tive hybrid functional we use PBEh 31 (PBE hybrid, also 
refer to as PBE1PBE 32 and PBEO 33 in the literature). 
For comparison, we also present results for the GUHF 
case. 

For Cr3 and Crs, we were able to verify that the two 
chiral magnetic structures (Fig. QJ, and Fig. [TJd, respec- 
tively) have the same total energy. These two chiral mag- 
netic states can be thought of as a product of a reflection 
of the spin density pseudovector in a molecular symmetry 
plane. For Cr 3 and Cr 5 , we have found that starting the 
SCF procedure from different initial guesses always leads 
to coplanar spin densities, although the plane containing 
the spin density does not necessarily coincide with the 
plane containing the nuclei since the spin density can ar- 
bitrarily rotate without changing the total energy. We 
therefore have chosen to constrain the spin magnetiza- 
tion to the plane containing the nuclei for the rest of our 
tests. 

In Figs. [2] and [3l we present a plot of the PBE spin 
density in the plane containing the nuclei for Cr 3 and 
Crs , respectively. The white (low spin polarization) holes 
at the nuclear positions are a consequence of the pseu- 
dopotential approximation. The red zones surrounding 
the nuclei correspond to high spin polarization regions. 
Four lobes can be distinguished around each atomic cen- 
ter, which is a signature of the spin polarization of the d 
orbitals. 

From Figs. [5]and [31 it can also be seen that the magne- 
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FIG. 2: Magnetization plot for Cr 3 obtained in a PBE calcu- 
lation. The arrows show the direction of the spin polarization 
(m/m) in the plane containing the nuclei, whereas the spin 
modulus m is represented in red. The top (a) and bottom (b) 
panels show two energetically degenerate configurations with 
+ and - chiralities, respectively. 



tization tends to be collinear in the atomic regions. Inside 
these atomic domains, the magnetization angle changes 
smoothly whereas it changes abruptly at the domain 
boundary. This was also observed in Fe clusters^ and in 
unsupported Cr monolayers in the 120° Neel state.— Spin 
density plots obtained using density functionals other 
than PBE do not differ qualitatively from these plots. 
However, the magnitude of the spin polarization does, as 
it is discussed below. 

In Table [II we summarize the results obtained for 
Cr clusters with the different functionals. In all cases, 
the atomic magnetization increases systematically when 
going from LSDA^PBE^TPSS^PBEh^GUHF. The 
value of (S 2 ) follows the same trend, and it can be taken 
as a measure of the total magnetization. As a remark, 
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(b) 

FIG. 3: Magnetization plot for Crs obtained in a PBE calcu- 
lation. The arrows show the direction of the spin polarization 
(m/m) in the plane containing the nuclei, whereas the spin 
modulus m is represented in red. The top (a) and bottom (b) 
panels show two energetically degenerate configurations with 
+ and - chiralities, respectively. 



we would like to recall that the Cr cluster geometries are 
fixed in these test calculations and hence relaxation ef- 
fects are not included in the reported atomic magnetic 
moments. For Cr3 and Cr5, we obtain comparable val- 
ues of the atomic magnetization m for a given functional. 
This is not the case for Cir and C112 clusters, where 
the atomic magnetic moments of the internal Cr atoms 
are smaller than those of the external atoms in all cases. 
From the values of m c and m e (see Fig. [T]) for Cry with 
PBEh and GUHF, we can notice a large effect of Hartree- 



FIG. 4: Magnetization plot for the Cr7 cluster obtained in 
a PBE calculation. The arrows show the direction of the 
spin polarization (m/m) in the plane containing the nuclei, 
whereas the spin modulus m is represented in red. 



Fock exchange compared to the rest of the clusters. In 
fact, GUHF gives an atomic magnetization in Cr7 ap- 
proximately 10 times larger than those obtained with 
LSDA, while for Cr3 and Crs the ratio is less than 2, 
and for &12 is about 3. 

In Figs. 2] and we present the PBE spin density in 
the plane containing the nuclei for Cry and Cri 2 , respec- 
tively. The plot for Cr7 resembles the spin density in 
the usupported Cr monolayer shown in Ref. Qjl The 
difference arises mainly in the low spin polarization re- 
gion of the Cr 7 cluster. Cri 2 can be thought as a cluster 
model of a two-dimensional Kagome lattice. However, as 
the magnetization of the internal and external atoms is 
not the same, one would expect that the actual ground 
state is a result of different competing effects. For in- 
stance, in a Heisenberg spin Hamiltonian, the uppermost 
Cr atom, Fig. [TJi, couples antiferromagnetically with its 
weaker polarized nearest neighbors and with its stronger 
polarized second nearest neighbors, and it is not clear a 
priori which coupling is larger. Therefore, we expect that 
this type of noncollinear DFT calculations would be help- 
ful to investigate the magnetic properties of clusters and 
molecules where a simple Heisenberg spin Hamiltonian 
cannot be straightforwardly applied. 



V. SUMMARY AND CONCLUSIONS 

We have generalized the treatment of the electronic 
spin degrees of freedom in density functional calculations 
to the case where the vector variables employed in the 
definition of the XC energy can vary in any direction. 
Our noncollinear generalization can be applied to gen- 
eral functionals containing a variety of ingredients. Our 
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generalization assumes that the XC energy depends on 
the local variables in the same manner as in the standard 
collinear case, and that the energy expression is invariant 
under rigid rotations of the spin quantization axis. This 
is not the most general way to define energy functionals 
for noncollinear magnetic systems, but it provides a gen- 
eral starting point to incorporate new terms like those 
suggested in Refs. andliHl. 

Test calculations on planar Cr clusters suggest that 
the choice of energy functional has an important impact 
on the resulting atomic magnetic moments, giving qual- 
itatively similar but quantitatively different results. We 
expect that our generalization will open the door to stud- 
ies on the performance of density functionals other than 
LSD A for noncollinear magnetic systems. 
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spin polarization (m/m) in the plane containing the nuclei, 
whereas the spin modulus m is represented in red. 



J.E.P thanks K. Capelle, R. Pino, A. fzmailov, and 
O. Hod for useful discussions and T. Van Voorhis for 
suggesting the Cri2 example. This work was supported 
by the Department of Energy Grant No. DE-FG02- 
0IERI5232, ARO-MUR1 DAAD-19-3-1-0169, and the 
Welch Foundation. 



1 U. von Barth and L. Hedin, J. Phys. C 5, 1629 (1972). 

2 P. Hohenberg and W. Kohn, Phys. Rev. B 136, 864 (1964). 

3 W. Kohn and L. J. Sham, Phys. Rev. A 140, 1133 (1965). 

4 Y. Tsunoda, J. Phys.: Condens. Matter 1, 10427 (1989). 

5 E. Sjostedt and L. Nordstrom, Phys. Rev. B 66, 014447 
(2002). 

6 D. Grohol, K. Matan, J.-H. Cho, S.-H. Lee, J. W. Lynn, 
D. G. Nocera, and Y. S. Lee, Nature Materials 4, 323 
(2005). 

7 G. Venturini, R. Welter, E. Ressouche, and B. Malaman, 
J. Alloys Compd. 210, 213 (1994). 

8 M. Uchida, Y. Onose, Y. Matsui, and Y. Tokura, Science 
311, 359 (2006). 

9 J. Kiibler, K.-H. Hock, J. Sticht, and A. R. Williams, J. 
Phys. F: Met. Phys. 18, 469 (1988). 

10 L. Nordsrom and D. J. Singh, Phys. Rev. Lett. 76, 4420 
(1996). 

11 T. Oda, A. Pasquarello, and R. Car, Phys. Rev. Lett. 80, 
3622 (1998). 

12 P. Kurtz, F. Forster, L. Nordsrom, G. Bihlmayer, and 
S. Blugel, Phys. Rev. B 69, 024415 (2004). 

13 S. Yamanaka, D. Yamaki, Y. Shigeta, H. Nagao, Y. Yosh- 
ioka, N. Suzuki, and K. Yamaguchi, Int. J. Quantum 
Chem. 80, 664 (2000). 

14 C. Kohl and G. F. Bertsch, Phys. Rev. B 60, 4205 (1999). 

15 R. C. Longo, E. G. Noya, and L. J. Gallego, Phys. Rev. B 
72, 174409 (2005). 

16 J. Mejfa-Lopez, A. H. Romero, M. E. Garcia, and J. L. 
Moran-Lopez, Phys. Rev. B 74, 140405(R) (2006). 

17 J. P. P. A. Ruzsinszky, J. Tao, V. N. Staroverov, G. E. 
Scuseria, and G. I. Csonka, J. Chem. Phys. 123, 062201 



(2005). 

18 A. D. Becke, J. Chem. Phys. 98, 5648 (1993). 

19 L. Kleinman, Phys. Rev. B 59, 3314 (1999). 

20 J. A. Pople, P. M. W. Gill, and B. G. Johnson, Chem. 
Phys. Lett. 199, 557 (1992). 

21 K. Capelle, G. Vignale, and B. L. Gyorffy, Phys. Rev. Lett. 
87, 206403 (2001). 

22 K. Knopfle, L. M. Sandratskii, and J. Kubler, Phys. Rev. 
B 62, 5564 (2000). 

23 R. McWeeny, Methods of Molecular Quantum Mechanics 
(Academic Press, San Diego, CA, 1989). 

24 Gaussian Development Version, Revision E.05, M. J. 
Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. 
Robb, J. R. Cheeseman, J. A. Montgomery, Jr., T. Vreven, 

G. Scalmani, K. N. Kudin, S. S. Iyengar, J. Tomasi, V. 
Barone, B. Mennucci, M. Cossi, N. Rega, G. A. Petersson, 

H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, 
J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Ki- 
tao, H. Nakai, X. Li, H. P. Hratchian, J. E. Peralta, A. 

F. Izmaylov, E. Brothers, V. Staroverov, R. Kobayashi, J. 
Normand, J. C. Burant, J. M. Millam, M. Klene, J. E. 
Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, 
R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, 
R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. 
Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. 

G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, 
O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, 
J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clif- 
ford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, 
P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, 
M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challa- 



9 



combe, W. Chen, M. W. Wong, and J. A. Pople, Gaussian, 
Inc., Wallingford CT, 2006. Gaussian, Inc., Pittsburgh PA, 
2003. 

P. Pulay, Chem. Phys. Lett. 73, 393 (1980). 

K. N. Kudin, G. E. Scuscria, and E. Cances, J. Chem. 

Phys. 116, 8255 (2002). 

M. Dolg, U. Wedig, H. Stoll, and H. Preuss, J. Chem. Phys. 
86, 866 (1987). 

S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys 58, 
1200 (1980). 

J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. 



Lett. 77, 3865 (1996). 

J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria, 
Phys. Rev. Lett. 91, 146401 (2003). 

J. P. Perdew, M. Ernzerhof, and K. Burke, J. Chem. Phys. 
105, 9982 (1997). 

M. Ernzerhof and G. E. Scuseria, J. Chem. Phys. 110, 
5029 (1999). 

C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 
(1999). 



